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ABSTRACT 


There  has  been  considerable  interest  lately  in  the  application  of  singular 
value  analysis  in  systems  theory.  The  basic  ideas,  however,  were  developed  in 
statistics  (Hotelling  introduced  principal  component  analysis  in  1933)  and  are 
currently  used  in  numerical  analysis  and  digital  filtering.  The  fundamental 
results  underlying  principal  component  analysis  are  presented  in  this  paper, 
and  these  results  are  applied  to  the  problem  of  tuning  multivariable  propor¬ 
tional  plus  integral  controllers.  Although  the  tuning  method  proposed  is  pre¬ 
liminary,  it  is  designed  to  avoid  possible  traps  which  would  prevent 
"tight"  tuning  with  conventional  tuning  of  individual  loops.  When  applied  to 
non-interacting  loops,  the  method  reduces  to  conventional  tuning  of  the  loops 
simultaneously. 
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I. 


INTRODUCTION 


Although  there  has  been  considerable  interest  lately  in  the  application 
of  singular  value  analysis  in  systems  theory,  [l]-[9],  the  basic  analysis 
techniques  involved  are  at  least  46  years  old.  Hotelling  [10],  [11],  intro¬ 
duced  principal  component  analysis  in  1933,  and  the  effectivenss  of  these  tech¬ 
niques  was  enhanced  substantially  by  the  development  of  an  algorithm  (S.V.D.), 

[12],  for  efficient,  accurate  computation  of  the  important  objects.  These 
techniques  are  currently  used  in  the  numerical  analysis  [13]-[1S],  and  digital 
filtering  [16]-[18].  Dempster  [19]  gives  an  excellent  geometric  treatment  of 
principal  component  analysis  as  well  as  an  overview  of  its  history  and  relation¬ 
ship  to  least  squares  approximation. 

In  this  paper  the  fundamental  results  underlying  these  analysis  tools  are 
presented  (Section  II),  and  a  preliminary  method  (using  these  tools)  for  on-line 
tuning  of  multivariable  proportional -plus -integral  controllers  is  proposed 
(Sections  III,  IV).  When  the  method  is  applied  to  a  set  of  non-interacting  loops, 
it  reduces  to  a  standard  classical  technique  ([20],  page  330). 

Although  the  point  will  not  be  pursued  here,  linearity  does  not  play  an 
essential  role  in  principal  component  analysis.  Preliminary  ideas  on  the  applica¬ 
tion  of  these  tools  to  nonlinear  problems  are  given  in  [21] . 


Notation:  IR  ,  will  represent  the  fields  of  real  and  complex  numbers.  For  a 

T  T  *  * 

vector  v  and  matrix  M,  v  ,M  represent  the  transpose,  v  ,M  represent  the  complex 
conjugate  and  v^,M^  represent  the  conjugate  transpose. 
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II.  FUNDAMENTAL  RESULTS 

The  results  given  in  this  section  form  the  foundation  for  the  principal 
component  analysis  techniques  to  be  discussed  in  later  sections.  The  two 
propositions  given  below  may  be  viewed  as  one  result  stated  in  two  different 
contexts:  the  first  involving  discrete  data  samples;  the  second,  continuous 
data.  The  constant  K  plays  no  role  other  than  an  aid  for  discussion. 

Discrete  Data: 

For  convenience  the  samples  are  organized  into  a  sequence  of  vectors 
y(l)  ,y(2) •  ,y(N)  in  <fn,  with 

N2  =  K  l  y(£]yH^.  (K>0) 

1=1 

This  matrix  is  positive  semidefinite  with  a  set  of  non-negative  eigenvalues 

2  2  2 
a.  z  a0  £  ...  2  0 

12  n 

and  corresponding  mutually  orthogonal  unit  eigenvectors  u1 ,u2> . . . ,un- 

Proposition  1A: 

Let  the  scalar  sequence  y^(l),y^(2), — ,y. (N)  be  defined  by 
Xitf)  -  ujy(£) 

n 

for  1  s  i  s  n.  Then  y(£)  =  l  u  y .(£),  where 

i=l  1  1 


r-*-* 


N 


k  I  |y±C^I 

£=i 


1/2 


=  o. 

1 


Proof;  The  first  identity  follows  from  the  fact  that  U=(u^ .u^, • . . ,un)  is  unitary. 
Furthermore 


K  I  llXi C-^D  j] 2  =K  I  u^y (•£) y  (£) Hu.  =  uVu. 

£=1  1  1=1  1  ill 

and  the  second  identity  follows  trivially.  □ 


Stated  in  other  words.  Proposition  1A  says  that  we  can  decompose  the  vector 
sequence  into  the  vector  sum  of  spacially  orthogonal  sequences  (called  components) 
ordered  with  respect  to  magnitude. 

The  important  objects  in  Proposition  1A  are  the  pairs  (c^.u^),  i=l,...,n, 

2 

and  these  can  be  computed  without  first  computing  W  .  If  a  data  matrix  Y  is 
constructed  with  N  columns  consisting  of  the  vectors  y(l) ,y (2) , . . . ,y(N) ,  then 


t 


and  (a^,u^),  lsisn,  are  the  singular  values  and  left  singular  vectors  of  Y.  These 
may  be  computed  directly  using  the  well  known  algorithm  (S.V.D.)  developed  by 
Golub  and  Reinsch  [12]  (see  also  [  9  ]). 


Continuous  Data: 


Now  consider  a  piecewise  continuous  map  f:$-*-Cn  defined  by  y=f(x),  and  let 

W2  -  K  f  2  fwAxJdx  (K>0) 

Jxl 

2  2  2 

with  eigenvalues  s  i  ^  on  i  0  and  mutually  orthogonal  unit  eigenvectors 

u.,u_,...,u  . 
i  i  n 

Proposition  IB: 

Let  be  defined  by  f^(x)=u^f  (x) .  Then 

f(x)=  l  U.f.(x) 
i=l  1 

where 

IC  f  2  |f  (x)|2dx 
X1 

Proof :  Similar  to  the  proof  of  Proposition  1A.  □ 

Here  we  decompose  the  vector  valued  function  into  the  sum  of  spacially  ortho¬ 
gonal  component  functions  ordered  with  respect  to  magnitude.  The  same  computational 

2 

tool,  S.V.D.,  can  be  used  to  compute  (a^ ,u^)  without  actually  computing  W  .  To 
do  this,  divide  [x^.Xj]  into  N  evenly  spaced  sample  points  and  construct  a  data 
matrix  Y  whose  columns  are  the  samples.  For  N  large,  rectangular  approximation 
of  the  integral  gives 


|  2  fCxJf^xJdx  a  (]/N)YYH 
,Jtl 

To  compute  C0^.^),  lsisn,  one  may  apply  S.V.D.  to  the  scaled  data  matrix 
(K/N) 1/2Y 


A  Combined  Result  for  Linear  Systems 

With  linear  systems,  we  shall  often  encounter  a  matrix  F(x)  which  can  be 
viewed  as  a  set  of  maps  ^ (x) ,f2(x) , . . . ,fm(x) ,  one  corresponding  to  each  column 
of  F.  In  this  situation  it  is  appropriate  to  combine  Propositions  IA,  IB.  Let 

w2  =  K  j  2  fi OO^iCxJdx  =  K  |  2  F(x)FH(x)dx. 

1_  X1  X1 

2  2 

If  (Oj.Uj) . . »(an,un)  are  the  ordered  eigenvalue,  eigenvector  (orthonormal) 
sets  for  W2,  and 

*i00  -  u“f(x) 


then 


where 


F(x)  -  I  uJf.(x) 
i=l  1  1 
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To  compute  (o^,u^)  in  this  situation,  divide  [x^,x2J  into  N  evenly  spaced 
sample  points  s^ ,s2> . . . .s^,  and  let 

Y  -  (F(sp  F(s2)  ...  F(sn)). 

If  there  are  m  columns  in  F  and  N  sample  points,  then  Y  has  mN  columns. 


Remarks  about  Perturbations: 

Suppose  the  matrix  F(x)  is  perturbed  to  give  FA(x)=F(x)+AF(x) .  Then 
u“fa(x)  =  Fi(x)  ♦  u“aF(x) 
and  it  follows  that 


’  ’ 
[X2||F.(x)-u“fa(x)  ||2dt 

1/2 

< 

f2  l|AF(x)  |[  2  dx 

*  X1 

1  J 

*  X, 

1  1 

where  the  equality  is  achieved  if  AF(x)  is  aligned  with  u^. 

This  is  a  double  edged  sword.  First,  it  gives  a  tool  for  coping  with  structural 
instability  associated  with  many  theoretical  results.  Theory  says,  for  example, 
that  if  every  column  of  F(x)  is  contained  in  a  proper  subspace  S  for  all  xe[Xj,x2], 
then  we  may  project  onto  S  to  simplify  the  situation  (by  reducing  dimensionality). 

Such  a  subspace  S  may,  however,  be  structurally  unstable:  there  may  exist  arbitrarily 
small  admissable  perturbations  such  that  the  columns  of  F(x)+AF(x)  are  not  contained 
in  a  proper  subspace.  We  are  guaranteed,  however,  that  F^(x)  will  have  weak  components 
(assuming  small  perturbations)  and  that  the  strong  components  will  define  a  subspace 
which  is  close  to  S. 


The  other  edge  of  the  sword  gives  us  help  in  deciding  roughly  how  accurate 
we  can  expect  the  components  to  be.  If  for  example,  there  are  components  whose 
magnitudes  are  of  the  same  order  as  the  instrumentation  precision,  then  they 
may  be  in  error  by  *<100%  and  one  can  hardly  use  them  (say  for  feedback)  with 
confidence. 


General  Comments  about  Applications: 

The  results  of  this  section  provide  a  strong  tool  for  spacial  analysis 
(possibly  on-line  in  many  situations)  of  multivariable  signals  in  the  time  or 
frequency  domains. 

If,  for  example,  f(t)  is  a  vector  of  signals  and 


N 


2 


f(t)fT(t)dt 


then  Proposition  IB  gives  a  decomposition  into  components  ordered  with  respect 
to  their  energy  (K=l)  or  average  power  K  = 


1 

f?(t)dt 


1/2 


a. 

i 


For  a  vector  f(jw)  in  the  frequency  domain,  and 


W2 


f(jw)fH(jw)dw, 


one  gets  information  about  the  spacial  distribution  of  f(jw)  over  the  frequency 


Remarks  about  Computation: 

Although  one  might  be  alarmed  at  the  thought  of  using  a  minicomputer  to 
compute  the  singular  values  and  left  singular  vectors  of  a  matrix  Y  with,  say, 

10  rows  and  100  columns,  it  is  quite  reasonable.  One  can  recursively  (treating 
one  column  of  data  at  a  time)  reduce  Y  to  a  unitarily  equivalent  matrix  (see  [14] , 
p.  383) 

A 

Y  =  YQ  (Q  unitary) 

A  _  2 

where  Y  is  10*10;  this  process  requires  «100(10  )=10,000  operations.  Singular 

a  3 

value  decomposition  of  Y  requires  k6(10  )=6000  operations,  giving  a  total  of 
operation  count  of  «16,000. 

At  first  it  may  seem  simpler  to  find  the  components  by  computing  the  eigen¬ 
values  and  eigenvectors  of  YY^.  Experts  recommend  that  this  be  avoided  (by 
using  S.V.D.)  for  the  following  reason  (see  [14],  p.  382  for  more  discussion): 

This  method  doubles  the  demand  for  internal  computer  resolution  associated  with 
the  algorithm. 

Specifically,  suppose  one  has  invested  money  in  12  bit  A-D  converters  and  has 
interfaced  these  properly  so  that  there  is  12  bit  resolution  associated  with  the 
data.  To  get  this  same  resolution  on  the  span  [0,a^] ;  that  is,  to  have 

-12  -12 
Oi  -  2  s  (computed  value  of  ck}  s  ck  +  2 

requires  at  least  12  bit  internal  resolution  using  S.V.D. ,  and  at  least  24  bit 
internal  resolution  using  the  "squared  up"  version  where  YY  is  computed. 


—  n  - 


III.  MULTIVARIABLE  PI  CONTROL 

Consider  the  following  PI  (Proportional* Integral)  control  loop  which  is 
assumed  to  be  open  loop  stable: 


With  a  single  isolated  loop  such  as  this,  one  can  often  follow  a  simple 
procedure  to  adjust  the  gains  k^,  k^;  i.e.  to  "tune"  the  control  loop.  One 
classical  method  (see  [20],  p.  330)  is  the  following: 

•with  kj=0  increase  k^  (this  improves  response  speed)  until  the  step  response 
is  highly  oscillatory.  Reduce  k^  by  a  factor  of  2. 

•increase  kj  (this  reduces  offset)  until  the  step  response  is  highly  oscillatory. 
Reduce  kj  by  a  factor  of  2. 

With  multiple  interacting  PI  control  loops,  there  are  inherent  "traps"  asso¬ 
ciated  with  extending  simple  procedures  such  as  this.  To  bring  some  of  these 
problems  into  closer  view,  consider  the  system  shown  below 


where  K^,  Kj  are  diagonal  matrices  of  proportional  and  integral  gains,  respectively. 
Let's  assume  that  each  variable  is  scaled  so  that  one  unit  corresponds  to  a  fixed 
percentage  of  "full  swing"  and  that  we  are  free  to  sample  e(t),  the  error  vector. 

Let  E(t)  be  the  matrix  made  up  of  error  responses  to  unit  steps;  i.e.  the 
il  column  of  E(t)  is  the  error  response  to  a  unit  step  applied  to  the  i  reference 
input.  For  a  linear  system,  the  response  to  an  arbitrary  vector  of  steps 


r (t)  =  r  6(t) 


is  given  by 

e(t)  =  E(t)r. 

The  steady  state  error  (assuming  stability)  is 

e  =  E  r 
ss  ss 

where  E  =  limit  E(t). 

ss  .  ' 

t-X» 

A  major  trap  follows  from  the  fact  that  an  operator  sees  only  proj ections 
of  the  vector  e(t)  on  the  basis  vectors  of  a  fixed  coordinate  system  associated 
with  the  physical  arrangement  of  hardware  (sensors) .  It  is  possible  for  rather 
simple  mechanisms  to  appear  complicated  in  this  fixed  coordinate  system.  The 
following  paragraphs  show  that  simple  mechanisms  involving  the  notions  of 
settling  time,  oscillations,  and  steady  state  errors  may  be  confusing  when 
viewed  through  projections. 


Settling  Time: 


Let  T(t)  =  E(t)  -  Egs  be  the  transient  response  map  so  that 


e(t )  =  etr(t) 


+  e 


ss 


where 


etr(t)  =  T (t)r 


It  is  certainly  possible  for  the  response  to  have  one  sluggish  component  which 
projects  significantly  onto  each  coordinate!  i.e.  all  loops  appear  sluggish. 

loop  2  error 


follows:  Let 


W2(t)  =  £  T(T)TTOOdT 


2  2  2 

with  eigenvalue,  eigenvector  pairs  (c^ (t) ,u1 (t) ) , (o2 (t) ,u2 (t)) , . . . , (o„ (t ) ,un(t)) 


t  .  =  minimum  time  such  that  a.  ft  .)  <  e. 
si  si' 

Note  that  t  ,  2  t  _  s  . . .  k  t  . 

si  s2  sn 

Further  on  in  this  section  we  shall  propose  one  procedure  for  "tightening" 
the  response  in  the  direction  associated  with  sluggish  components. 

Steady  State  Errors: 

2  A  T  2 

Consider  W  =  E  E  with  eigenvalue,  eigenvector  pairs  (a  . ,u  ),..., 
ss  SS  SS  SS i  S S J. 

2 

(0ssn»ussn^‘  **  *s  entirely  possible  that  there  are  only  a  few  strong  offset 
error  components  which  project  onto  every  loop 


Oscillations : 


To  avoid  unnecessary  complications,  let's  assume  that  there  is  a  single 

lightly  damped  sinusoidal  component  observed  in  the  loop  error  responses,  and 

f2ir/w 


that  the  observed  frequency  is  w.  Let  T 


av 


T (t)dt  and 


2it/w 

0 


(T(t)-Tav) (T(t)-Tav)‘dt 


2  2 

with  eigenvalue,  eigenvector  pairs  (°trl >utrl) » • • • » Cotrn»utrn) 


-iff 


Again  it  is  possible  that  there  is  one  strong  oscillatory  component  which 


projects  onto  a  number  of  loops 


loop  2 


-a  j,u  j  --  corresponds  to  the  principal 
oscillatory  component 


loop  1 


A  tuning  method  which  allows  one  to  avoid  these  traps  is  given  in  the  next 
section. 


IV.  A  PRELIMINARY  ON-LINE  TUNING  METHOD 

These  ideas  are  preliminary  and  essentially  untested.  Undoubtedly  tests 
currently  in  progress  will  lead  to  better  tuning  methods  than  the  rather  vague 
one  proposed  in  the  following  paragraphs.  The  tuning  method  is  basically  the 
classical  technique  (given  in  Section  III)  applied  to  principal  components.  If 
loops  are  noninteracting,  the  method  reduces  to  the  classical  technique  applied 
to  all  loops  simultaneously. 

Part  1  Tightening  the  Response: 

Let's  assume  that  initially  Kj=0,  and  that  the  diagonal  components  of  Kp  are 
increased  (at  the  same  rate,  if  you  like)  until  there  is  a  highly  oscillatory 
component.  If  the  components  of 


T(t)  -  T 


av 


have  magnitudes  atrj > • • • >atrn  are  roughly  equal,  then  this  step  is  essentially 

complete;  reduce  the  gains  by  a  factor  of  two  (to  reduce  the  oscillations),  and 
proceed  to  Part  2. 

If  there  are  some  weak  components,  then  increase  gains  in  these  directions. 
Specifically,  suppose  »  °trk+1 >  and  let 


U1  '  (U1  u2 


"k); 


U2  ’  <Vl 


Then  in  the  proportional  branch,  we  insert 


where  a,B  are  tuning  parameters  with  a=B=l  initially.  The  parameter  a  should  be 
increased  until  oscillations  appear  in  the  lower  branch;  it  may  be  necessary  to 
reduce  g  in  the  process  to  maintain  stability;  i.e.  to  keep  upper  branch  oscillations 
bounded. 

This  process 

•compute  components  of  T(t)-T 

•insert  a  Rotation-Tuning  Gains -Rotation  block 

•adjust  gains  to  get  oscillations  in  weak  branch 

can  be  continued  until  the  oscillations  are  "full"  with  no  weak  components.  The 

A 

resulting  gain  matrix  Kp  in  the  proportional  branch  is  the  product  of  Kp  and  the 

A 

inserted  blocks.  The  elements  of  Kp  should  be  reduced  by  a  factor  of  2  to  reduce 
oscillations.  This  completes  Part  1. 


Part  2  --  Inserting  integral  gain  to  reduce  offset  errors: 

At  this  point  it  may  be  necessary  to  insert  integral  action  if  there  is 

A 

significant  offset.  Here  we  shall  use  the  components  of  Egs  (with  controller  K) 

If  o  .  »  o  ,  , ,  let 
ssfk  ssk+1 


^1  ^Ussl  Uss2  Ussk^ ’  ^2  ^Ussk+1  Ussn^ 


and  configure  the  controller  in  the  following  way: 


with  this  configuration,  a  should  be  increased  from  zero  until  oscillations  occur 
If  the  oscillations  are  "full"  in  the  space  corresponding  to  the  integral  branch, 
reduce  a  by  a  factor  of  two  and  stop. 

If  the  oscillations  in  the  integral  branch  are  not  full,  we  may  proceed  as 

follows  iteratively  until  the  oscillations  are  full. 

T 

•compute  components  of  (T(t) -Tav)Uj 

•insert  a  Rotation-Tuning  Gains-Rotation  block  (as  in  part  1)  in  the  integral 
branch 

•adjust  gains  to  increase  oscillation  in  weak  integral  branch 

We  are  one  step  away  from  completion  of  the  proposed  tuning  algorithm,  with 


the  following  controller  structure 


The  last  step  is  to  reduce  the  elements  of  Kj  by  a  factor  of  two. 


Concluding  Remarks: 

It  is  certainly  possible  in  many  process  control  problems  to  use  a  single 
microprocessor  system  (<$10,000)  for  multivariable  PI  control  of  a  small  number 
of  loops  (say  10) .  Process  control  systems  of  this  type  have  been  applied 
successfully  (see  [22],  [23])  in  the  control  of  industrial  heating  and  air 
conditioning  equipment  (boilers,  chillers,  cooling  towers,  air  handlers,  etc.). 
The  tuning  method  proposed  in  this  section  can  be  implemented  with  little  or  no 
additional  hardware.  The  procedure  is  simple  and  consistent  with  classical 
tuning  methods..  With  a  well  designed  man-machine  interface  (say  bar  graphs  for 
singular  values  and  automatic  generation  of  rotation  blocks)  the  method  would 
probably  be  acceptable  to  plant  operators. 
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ROBUST  STABILITY  OF  LINEAR  SYSTEMS  -  SOME 
COMPUTATIONAL  CONSIDERATIONS* 


by 


Alan  J.  Laub** 


1.  INTRODUCTION 

In  this  paper  we  shall  concentrate  on  some  of  the  computational  issues 
which  arise  in  studying  the  robust  stability  of  linear  systems.  Insofar 
as  possible,  we  shall  use  notation  consistent  with  Stein's  paper  [1]  and 
we  shall  make  frequent  reference  to  that  work. 

As  we  saw  in  [1]  a  basic  stability  question  for  a  linear  time-invariant 
system  with  transfer  matrix  G(s)  is  the  following:  given  that  a  nominal 
closed- loop  feedback  system  is  stable,  does  the  feedback  system  remain 
stable  when  subjected  to  perturbations  and  how  large  can  those  perturba¬ 
tions  be?  It  turned  out,  through  invocation  of  the  Nyquist  Criterion, 
that  the  size  of  the  allowable  perturbations  was  related  to  the  "nearness 
to  singularity"  of  the  return  difference  matrix  I  +  G(jto) .  Closed- loop 
stability  was  said  to  be  "robust"  if  G  could  tolerate  considerable 
perturbation  before  I  +  G  became  singular. 
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We  shall  now  indulge  in  a  modicum  of  abstraction  and  attempt  to 
formalize  the  notion  of  robustness.  The  definition  will  employ  some 
jargon  from  algebraic  geometry  and  will  be  applicable  to  a  variety  of 
situations.  While  no  deep  results  from  algebraic  geometry  need  be  em¬ 
ployed,  the  exercise  of  formulating  a  precise  definition  is  a  useful  one 

for  clarifying  one's  thinking. 

N 

Let  p  €  3R  be  a  vector  of  parameters  from  some  problem  being  studied 
and  suppose  we  are  interested  in  some  property  II  of  this  data.  The  vector 
p  may  consist  of  the  elements  of  various  matrices,  for  example.  If  II 
is  true  at  some  nominal  parameter  set  pQ  we  are  frequently  concerned  with 
whether  II  remains  true  in  a  "neighborhood"  of  pQ . 

For  example,  pQ  may  be  the  elements  (a^,  ...,  a^,  a^,...,  &nn 
of  a  nonsingular  nxn  matrix  AQ  and  we  are  interested  in  the  nonsingularity 
of  nearby  matrices.  We  shall  proceed  to  formalize  the  often-heard  statement 
that  "almost  all  nxn  matrices  are  nonsingular".  First,  the  jargon: 

Definition  1:  A  variety  1/  *  {p  e  1RN :  4^  (p^, . .  • ,  PN)  =  0,  i  =  1,...,  k} 

where  (x, , . . . ,  x„)  e  3R[x, , . . . ,  x  1  are  polynomials . 

1  1  N  1  N 

N 

\J  is  proper  if  !/  and  nontrivial  if  (/  f  <p. 

N  r  1 

Definition  2:  A  property  is  a  function  II:  3R  -*•  {0,  1).  The  property 
TI  holds  if  IT  (pi  =  1  and  fails  if  H(p)  =0, 

Definition  3:  If  V  is  a  proper  variety,  fl  is  generic  relative  to  1/ 
provided  II  (p)  =0  only  if  p  e  I/.  A  property  H  is 
generic  if  such  a  !/  exists. 

Our  discussion  to  this  point  is  purely  algebraic.  Now  let  us  intro¬ 
duce  a  topology  onlRN,  say  the  topology  induced  by  some  vector  norm  jj  •  ||  . 
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Furthennore,  let  1/  be  any  nontrivial,  proper  variety.  Then  we  have 
the  following  topological  definition. 

Definition  4:  The  property  II  is  well-posed  at  p  e  t/C(the  complement  of 
V)  if  II  also  holds  in  a  sufficiently  small  neighborhood 
of  p. 

Lemma  1;  The  set  S  of  points  where  a  generic  property  is  well-posed 

•  c 

is  open  and  dense.  Moreover,  the  Lebesgue  measure  of  S 

is  zero. 

The  proof  of  Lemma  1  is  routine  and  is  omitted.  It  is  easy  to  see 
that  a  point  p  where  a  generic  property  holds  is  well-posed  but  that  the 
converse  is  not  necessarily  true. 

We  now  have  sufficient  framework  to  make  a  formal  definition  of 
robustness. 

Definition  5:  Given  a  point  p  with  generic  property  II  (generic  with 

respect  to  some  proper  variety  V)  well-posed  at  p,  let 

d  =  min  [|  p  —  v  ||  . 
vet/ 

We  say  II  is  robust  at  p  if  d  is  "large". 

The  number  d  is  frequently  difficult  to  compute  or  estimate.  When 
it  can  be  determined,  it  gives  valuable  information  about  how  much 
perturbation  or  uncertainty  can  be  tolerated  at  p.  For  the  situation 
of  special  interest  in  this  paper, Example  2  below,  we  shall  see  that 
d  can  be  explicitly  calculated,  at  least  theoretically.  We  now  illustrate 
the  above  concepts  with  two  examples. 


This  example  is  chosen  from  Wonham  [2]  who  uses  the  concepts  of 


genericity  and  well-posedness  in  nontrivial  ways  for  a  variety  of  control- 
theoretic  problems.  In  this  trivial  example,  we  seek  solutions  of  the 
system  of  linear  equations 


Ax  *  b 

where  A  e  IR0301  (i.e.,  A  is  an  mxn  matrix  with  real  coefficients)  and  b  e  IR™. 

Our  parameter  vector  is  p  where 

T  N 

p  «  (a  /...,  a  ,  a  ;  b  ,...,  b  )  6  m  ,  N  =  mn  +  ra 

11  In  mn  1  m 


(  denotes  transpose) .  II  is  the  property  of  the  equation  having  a  solution 

which  is  equivalent,  of  course,  to  the  statements  that  b  6  Im  A  or 

(l  2\  jb\ 

rk[A,  b]  =  rk  A.  For  example,  if  A  =  I  and  b  =1  I  then 


(0  if  b  *  2b. 
H(l, 2, 2, 4;  b  -h  )  =  ) 


1'  2 


(  1  if  b,  -  2bl 


It  is  then  easy  to  show  the  following:  (see  [2]) 


1.  II  is  generic  if  and  only  if  m  <_  n. 

2.  II  is  well-posed  at  p  if  and  only  if  rk  A  =  m. 


Example 


This  example  is  similar  to  Example  1  in  the  special  case  ra  =  n.  We 


are  given  a  nonsingular  matrix  A  S  3RnXn  and  we  are  concerned  with  the 


nearness  of  A  to  singularity.  Identifying  A  with  p  =  (a^»' 

a_,  , ...»  a  )  we  define  the  property  II  by 
21  nn 


in' 
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!0  if  p  represents  a  singular  matrix 

1  if  p  represents  a  nonsingular  matrix  . 

Then  it  is  easy  to  see  that  II  is  a  generic  property  and  well-posed  where 
it  holds.  This  is  the  precise  statement  that  "almost  all  nxn  matrices 
are  nonsingular".  Formally  writing  down  the  determinant  of  A  as  a  poly¬ 
nomial  in  a, a  defines  the  necessary  variety  V .  It  turns  out, 

11  nn 

in  a  theorem  attributed  by  Kahan  [3]  to  Gastinel,  that  the  distance  d 

( 2 

from  a  point  p  6  V  to  1/  can  be  explicitly  determined. 

Theorem  1:  A  nonsingular  matrix  A  differs  from  a  singular  matrix  by  no 
more  in  norm  than  — — ,  i.e.,  given  A, 

II*-1  II 


=  min{  ||  E  |1  j 


A  +  E  is  singular}  . 


Thus  d  *  - ; —  and  we  might  say  that  A  is  robust  with  respect  to 

II*'1  II 

invertibility  if  d  is  "large".  To  avoid  certain  scaling  difficulties, 
it  may  be  more  desirable  to  work  with  a  relative  measure  of  distance, 
drel,  defined  by 


^rel  _  d  _  _ 1 _  _  1 

'  11*11  ‘  11*11  •  ll*"1!!  "  <w 

The  quantity  k(A)  is  recognizable  as  the  condition  number  of  A  with 
respect  to  inversion.  Of  course,  all  the  above  quantities  depend  on  the 
particular  matrix  norm  used.  To  exhibit  the  specific  dependence  on  the 
norm  ||  •  ||  we  shall  append  a  subscript  "q".  For  example, 
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*  --T-  • 

q  II*"1  n. 

The  minimizing  E  in  Theorem  1  can  be  explicitly  constructed  for  a  number 
of  standard  matrix  norms.  For  example: 

1-  H*ll2  ■  <W*T*,,l/2  • 

T  ,  „  rixn 

Let  A  have  singular  value  decomposition  A  =  USV  where  U,  V  SIR 

are  orthogonal  and  S  =  diag{a, 0  }.  The  a.'s,  a  >...>  a  >  0, 

X  n  i  x  —  —  n 

T 

are  the  singular  values  of  A.  The  minimizing  E  is  given  by  E  *=  URV 
where  R=  diag{0,...,  0,  -o  }.  Then 


IIa"1!! 


and  A  +  E  is  singular.  The  singular  direction,  i.e.,  a  nonzero 
vector  z  such  that  (A  +  E)  z  =  0 ,  is  given  by  the  n—  column  of 
V. 


2.  ||a|| 


n 

““  (  l 

ien  j«l 


{l,2, ... ,  n)  . 


Suppose  A 


-1 


and 


-li 


11 

»  £  | a  |  for  ken.  Then  the 

j-1  D 


minimizing  E  is  a  matrix  all  of  whose  elements  are  0  except 


for  the  k—  column  which  consists  of  the  elements 


with  the  only 


0 
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th  T 

nonzero  component  of  u  being  in  the  k —  row,  we  have  E  »  -zu 


and  clearly  ||e||  = 


IK 


-li 


Now,  (I  +  EA  Sz  *  (1  -  UTA  *z)z“  0 


since  the  k—  element  of  A~^z  is  £  l^jl  =  IK  ^11 « 90  that  u  A  z  =»  1 . 

-1 

Hence  A  +  E  =  (I  +  EA  )A  is  singular.  Moreover,  the  singular  direction  is 
given  by  A  ^z  since  (A+E) A  ^z  =  0. 


T  -1 


n 


3. 


max 


^  '  il 

j6n  1=1  J 


U  kj}. 


and  can  be  derived 


The  results  for  this  norm  are  analogous  to  ||  • 

directly  or  *  noticic,  that  II  A  ||,  -  I  »’lL  .  For  co^leteaes,  « 

note  that  if  ||  A  ||^  =  £  laijJ  for  k  e  n  and 

i=l 


sgn 


u  = 


i/IK”1! 


then  the  minimizing  E  is  given  by  E  =  -uz  . 


We  shall  see  in  Section  3  how  the  results  in  Example  2  can  be  applied 
in  studying  robustness  of  stability  of  linear  systems. 
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2.  THE  LINEAR  SYSTEMS  SETTING 

In  this  section  we  shall  provide  a  brief  introduction  to  both  the 
linear  time-invariant  systems  setting  and  to  the  fundamental  notion 
of  feedback.  This  will  serve  a  two- fold  purpose:  first,  to  set  the  stage 
for  the  basic  stability  results  and  second,  to  introduce  the  jargon  and 
notation,  especially  for  non-engineers.  This  material  is  standard  and 
can  be  found  in  any  of  a  number  of  standard  textbooks  on  control  systems. 

We  shall  consider  modelling  physical  systems  by  models  which  take 
the  form  of  a  system  of  linear  constant-coefficient  ordinary  differential 
equations 


X (t)  «  Ax(t)  +  Bu(t) 


(1) 


y(t)  -Cx(t) 


(2) 


The  vector  x  is  an  n-vector  of  states ,  u  is  an  m-vector  of  inputs  or 
controls,  and  y  is  an  r-vector  of  outputs  or  observed  variables. 

Starting  from  the  initial  condition  x(0)  the  solution  of  (1)  is 
well-known  to  be 


0 /■ 


x(t)  -  etAx(0)  +  /e(t-T)A Bu(T)dT,  t  >  0 


(3) 


so  that  the  output  is  given  by 


y(t)  »  CetAx(0)  +  J ce^  T^ABu(T)dx,  t  >_  0 


(4) 


tA 


where  e  is  the  matrix  exponential  defined,  but  not  generally 
computed,  by 


tA 


+°°  k  k 

I  tA 


k-0 


kl 
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The  matrix  Ce  b  is  called  the  impulse  response  matrix. 

Denoting  (one-sided)  Laplace  trams forms  by  upper  case  letters,  take 
Laplace  transforms  in  (4)  to  get 

Y(s)  =  CX(s)  =  C(?I  -  A)-1x(0)  +  C(sl  -  A)-1BU(s)  .  (5) 

The  matrix  G(s):  =  C(sl  -  A)  ^"B  is  called  the  transfer  matrix.  Notice 
that  G(s)  is  the  Laplace  transform  of  the  impulse  response  matrix. 

As  will  be  seen  in  the  sequel,  it  is  of  interest  to  3tudy  the 
response  of  the  above  linear  system  to  sinusoidal  inputs  of  the  form 

u(t)  =■  e-^v,  t  o  (6) 

where  v  is  a  constant  m-vector,  (d  is  the  frequency  of  the  sinusoidal 
input,  and  j  =  /-T.  The  response  of  (1)  to  this  input  cam  then  be 
shown  to  be  of  the  form 

x (t)  »  etAa  +  {jcol- Aj^Bve^,  t  >_  0  (7) 

where  a  is  a  constamt  n-vector  depending  on  initial  conditions.  Now, 

in  the  case  where  A  is  staible  (i.e.,  its  spectrum  lies  in  the  left- 

tA 

half  of  the  complex  plane)  the  quantity  e  a  goes  to  zero  as  t  approaches 
+°°.  The  resulting  output 

y(t)  =  C(jajl  -  A)  1Bve^CUt  (8) 

is  called  the  steady-state  frequency  response  amd  the  matrix 

G ( jOJ)  :  *  C( jcol  -  A)"1B  ,  (9) 

which  turns  out  to  be  the  transfer  function  evaluated  at  s  *  jo>,  is 
called  the  frequency  response  matrix. 
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Turning  now  to  the  case  of  a  real  signal  given  by 

u^t)  =  v^sin  (cot  +  1>k) ,  t  ^  0  (10) 

u. (t)  =  0 ,  i*l,...,m;ij*k, 

we  have  steady-state  frequency  response  of  the  £th  output  given  by 


Yo  (t)  “  iGov(3“)|v.  sin  (o)t  +  <J>.  +  iK.) 


(11) 


-  1  Vjw/ 1  vk  'r  ^ 

where  i/>£k  -  arg(G)lk  (jcu) )  . 

Aside  from  its  obvious  importance  in  the  above  analysis,  the 
frequency  response  matrix  is  important  for  two  reasons: 

1.  Sinusoidal  signals  are  readily  available  as  test  signals 

for  a  linear  system  so  G(jai)  can  be  experimentally  determined. 

m 

2.  Various  plots  or  graphs  associated  with  G(jco)  can  be  used  to 
analyze  control  systems,  for  example,  with  respect  to  stability. 
Plots  such  as  those  associated  with  the  names  of  Bode,  Nichols, 
and  Nyquist  are  essentially  different  ways  of  graphically 
representing  lG£k(jw)|  and  arg(G^k  ( jcu) )  as  functions  of 

<o.  These  plots  are  used  extensively  in  the  analysis  of 
single-input  single-output  control  systems  where  the  robust¬ 
ness  of  stability,  e.g.,  the  amount  of  gain  and  phase  margin 
available,  is  checked  essentially  visually.  The  appropriate 
techniques  in  the  multiple-input  multiple-output  case  are 
still  being  investigated  and  part  of  the  motivation  for  the 
research  in  [1]  and  this  paper  is  directed  towards  this  end. 
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Tuming  now  to  the  notion  of  feedback  whose  essential  idea  is  to 
allow  for  stability  of  a  system  in  the  face  of  uncertainty  (noise, 
model  error,  etc.),  the  diagram  below  illustrates  the  basic  (unity) 
feedback  control  system: 


Fig.  1.  Basic  Feedback  Control  System 

Here  u  is  a  reference  input,  y  is  the  output,  and  e  =  u  -  y  is  the  error 
or  difference  between  the  reference  input  and  the  output  which  we  wish 
to  be, ideally,  zero.  The  plant,  compensators,  actuators,  and  sensors 
are  all  represented  by  G.  There  are  much  more  elaborate  and  detailed 
feedback  structures  than  that  described  above  and  the  structure  can  be 
studied  in  a  considerably  more  general  function-space  setting  (see  (4] , 
for  example)  than  the  simple  linear  causal  time-invariant  setting  we 
shall  consider.  However,  the  simple  system  is  adequate  to  exhibit  most 
of  the  key  ideas  in  this  paper.  Now,  in  this  system  we  have 

e  =  u-  y  =  u-  Ge  (12) 

or. 


(1  +  G)e  =  u 


(13) 
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The  quantity  I  +  G  is  called  the  return  difference  matrix.  As  in  [1] , 
the  matrix  G ( joi)  then  provides  sufficient  data,  via  the  Nyquist  criterion 
to  test  for  stability  of  the  closed-loop  system.  Henceforth,  we  shall 
assume  that  our  nominal  feedback  system  above  is  stable  in  which  case 
I+G  is  invertible.  Then  from  (13)  we  have 

e  «  (I  +  G)*”^u  (14) 


so  that 

y  =  Ge  »  G (I  +  G)  .  (15) 

In  (15),  the  quantity  G(s) (I  +  G(s))  ^  is  called  the  closed-loop  transfer 
matrix  while  G(jco)  (I  +  G(j&i))  is  called  the  closed-loop  frequency 
response  matrix.  We  then  pose  the  basic  stability  questions 


Does  the  nominal  feedback  system  remain  stable  when  subjected 
to  perturbations  and  how  large  can  those  perturbations  be? 


Let  us  observe  at  this  point  that  there  is  nothing  sacred  about 
linearity  in  the  above  discussion  and  more  general  nonlinear  treat¬ 
ments  can  be  found  in  [4]  and  [5],  for  example.  The  question  of  "near¬ 
ness  to  singularity"  of  (I  +  G),  even  in  the  nonlinear  case,  is  naturally 
intimately  related  to  a  notion  of  condition  number  for  nonlinear 
equations.  The  interested  reader  could  readily  adapt  the  ideas  of 
Hheinboldt  [6]  to  the  particular  application  at  hand  here. 
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3.  BASIC  STABILITY  RESULTS  AND  RELATED  TOPICS 

a.  ADDITIVE  AND  MULTIPLICATIVE  PERTURBATIONS 

We  shall  consider  two  fundamental  types  of  perturbations  in  the 
basic  feedback  system  of  Fig.  1.  Throughout  this  section,  [|  •  ||  will 
denote  any  matrix  norm  with  ||  I  ||  =  1.  The  first  case  to  be  considered 
is  the  case  of  additive  perturbations  to  G,  pictured  below: 


Fig.  2.  Additive  Perturbations 

In  other  words,  the  nominal  G  is  perturbed  to  G  +  L.  Under  the  assumptions 
that  both  the  nominal  closed-loop  system  and  the  perturbation  L  are 
stable  it  can  be  seen  from  the  Nyquist  criterion  and  the  identity 

I  +  G  +  L  =  (I  +  G)  [I  +  (I  +  G)-1LJ  (16) 

that  the  perturbed  closed-loop  system  remains  stable  if 

|(<I  +  G(ju>>r1L(ju>)  ||  <1,  «  >  0  (17) 

A  weaker  condition  them  (17)  but  one  which  directly  exposes  L  is 


I 


|  (I  +  G  ( ju) ) 


CD  >  0 


The  second  case  to  be  considered  is  that  of  multiplicative  perturba¬ 
tions: 


u  +  _  e 


1+  L 

G 

_ : 

Fig.  3.  Multiplicative  Perturbations 

In  this  case,  the  nominal  G  is  perturbed  to  G{I  +  L).  Under  the  assumptions 
that  both  the  nominal  closed- loop  system  and  the  perturbation  I>  are 
stable  it  can  be  shown  from  the  Nyquist  criterion  and  the  identity 


I  +  G(I+L)  =  (I  +  G)  [I  +  (I+G-1)_1L) 


that  the  perturbed  closed-loop  system  remains  stable  if 


(I  +  G_1(j<D))"1L(j«)  ||  <  1, 


w  >  0 


(assuming  G  exists).  Again,  a  weaker  condition  than  (20)  but  one 
which  directly  exposes  L  is 


L  (jw)  <  - = - - -  ,  0)  >  0  . 

||(I  +  G'1{ja)))“1|| 


Remark  1:  As  we  noted  in  Section  1,  the  above  inequalities  are  tight 


i.e.,  the  <  cannot  be  replaced  with  . 

Remark  2:  Where  convenient  we  shall  henceforth  drop  the  "jcu"  arguments. 
Remark  3:  It  must  be  stressed  that  the  results  based  on 

11(1+  G*1)'1!!  Hi.  II  <  1  (18).  (21) 

are  weaker  than  those  based  on 

||(I+  G±X)"lL  I!  <  1  (17) ,  (20) 

since 

||(I+G±1)“1l||  <  ild+G41)"1!!  •  II L  ||  .  (22) 

±1  ,  , 

For  example,  if  L  =  c(I+G  )  for  some  constant  c,  |c|  <  1,  the 
differences  in  the  bounds  are  obvious.  In  (18),  (21)  we  have 

||(I  +  G±1)"1||  •  ||  L  ||  =  |  C  |  •  K(I  +  G41) 
while  in  (17) ,  (20)  we  have 

||<t.5iV1L||  -  1  cl 
and  it  is  possible  to  have 

| C |  «  |c I  •  K  (I +  G-1) 

However,  for  random  perturbations  L,  (22)  is  often  approximately  an 
equality.  To  see  this,  note  that  a  random  (dense)  L  will  almost  surely 
be  invertible;  recall  Example  2.  It  is  then  easy  to  show  that 

||d  + g±1)_1l||  £  ||(I  +  G±1)"1||  •  ||  I.  ||  . 

\W  I 


Again,  since  L  is  random,  it  will  almost  surely  be  well-conditioned 
(w.r.t.  inversion)  so  that  [|l  •  Hence, 


11(1  + 


as 


A  related  aspect,  also  worth  noting,  follows  from  the  inequalities 

L  II  .  it .  _il.  “1.  ii  .  n  _il.  -In  ii  _  i 


(I +  G±X) _1 I 


+1 

K  (I  +  G  ) 


<  ||(I  +  G  )  l||  <  |j(I  +  G-  ) ' 


±1  ±1 
If  (I  +  G  )  is  reasonably  well-conditioned  (<(I+G  )  near  1)  ,  the 


majorization  (22)  will  not  be  a  bad  overestimate. 


Remark  4;  By  our  discussion  in  Section  1,  the  appropriate  measure  of 
stability  robustness  is 


d  *  min 

OJ>0 


_ 1 _ 

±1  -111  * 

||(I+G  X(joj) )  || 


(23) 


and  in  the  sequel  we  shall  consider  methods  of  efficiently  plotting 


as  a  function  of  oi. 


This  quantity  is  familiar  from 


classical  sensitivity  analysis  where  it  is  shown,  in  the  single-input 
single-output  case,  that  the  change  in  the  output  of  a  closed-loop 
system,  due  to  (additive)  perturbations  in  G  (scalar) ,  is  reduced  by 
a  factor  of  1  +  G  compared  with  the  open-loop  effect. 


Remark  5:  So  far  we  have  required  nothing  of  our  norm  other  than 
||  I  ||  ■  1.  Of  course,  a  frequently  occurring  norm  in  much  of  the 
analysis  of  linear  systems  is  the  spectral  norm  ||*  |L  .  In  that  case 
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||(I  +  G±1)‘1| 


±1 

is  the  smallest  singular  value  of  (I  +  G  ).  Let 


d  (04) 

q 


(I  +  G^jto))"1^ 


(24) 


Vfe  are  interested  in  plotting  d^(oj)  versus  a>  for  large  numbers  of 
aj's.  We  shall  see  in  the  sequel  that  determining  d^  (to)  can  be 
somewhat  more  expensive  to  determine  than,  say  d^  (co)  or  d^Cfo).  More¬ 
over,  note  that 


-T-IIMIjilUlljl-Ma 

vm 


(25) 


and 


~  II  A  112  i  II  A  IL  1  ^  II  A  |l2  (26) 

for  A  6  <rnutm.  Since  we  are  usually  most  interested  in  order-of-magnitude 

estimates  of  d  (co) ,  d_  (oj)  will  lie  in  a  strip  sufficiently  close  to 
q  2 

dx (04)  ,  for  example,  to  give  the  same  qualitative  information.  The 
number  m  which  is  the  number  of  inputs/outputs  in  the  system  is  typically 
no  more  than  about  10  and  is  frequently  much  less. 


b.  PFT.ATIONSHIPS  between  additive  and  multiplicative  perturbations 

The  following  theorem  relates  additive  and  multiplicative  perturba¬ 
tions.  Again,  the  "joj's"  will  be  omitted  for  convenience  and  all 
relations  will  be  assumed  to  hold  for  all  04  >  0. 
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Theorem  2: 


|(i  +  g'1)"1| 


(1+  G) 


-li 


<  1 


Proof ;  From  the  identity 


(I  +  G-1)"1  +  (I+G)"1  =  I 


(27) 


we  have 


Id  +  G-1)-1! 


|(I+  G) -1 1[  <  ||(I+  G-1)"1  +  (I  +  G)‘1||  =  III 


We  now  get  immediately  the  following  useful  corollary : 

Corollary  1:  Assuming  that  both  the  nominal  closed-loop  feedback 
system  of  Fig.  1  and  the  perturbation  L  are  stable  then  the  perturbed 
system  is  stable  under: 


(a)  additive  perturbations  if 
1 


L  < 


1  +  Ild+G'1)"1! 


(28) 


(b)  multiplicative  perturbations  if 


L  < 


1  +  (I  +  G) 


-li 


(29) 


Proof;  Follows  immediately  from  Theorem  2  noting  that 


i  +  I!  (i  +  g*1)"1  ||  Ild  +  G71)"1!! 


C.  SPECIAL  RESULTS  FOR  THE  SPECTRAL  NORM 

In  this  subsection  we  shall  present  some  results  related  to  those 
in  subsections  a.  and  b.  but  specialized  to  the  ||*  ||  - 


norm.  For 


L9- 


a  matrix  H  e  <E”-  with  singular  values  a  (H)  >  ...>  a  (H)  >  0  we  note 

1  —  —  m  — 

that  ||  H  j|2  =  cr^H).  If  H  is  nonsingular,  ||h  1  (| 2  =  >0.  In  the 

m 


•  -  norm  (28)  becomes 


0  (I  +  G_1) 

a  (L)  <  — - rr 

1  +  O  (I+G  L) 
m 


while  (29)  becomes 

a  (I  +  G) 

ai(L)  <  1  +  a  (I  +  G)  * 
m 

We  shall  make  great  use  in  the  sequel  of  the  following  result 
of  Fan  [7] . 


Theorem  3:  Let  A,B  6  <C  .  Then 


(a)  °i+j_1<A+B>  +  ^(B)?  i  >  1,  j  >_  1 


(b)  Oi+j_1(AB)  _<  ^(AjOj  (B);  i  >  1,  j  >  1 


Part  (b)  of  Theorem  3  can  be  used  to  relate  0  (I+G)  and  a  (I+G  ) 

m  m 

Theorem  4:  (a)  — V — 0  (I  +  g"1)  <  a  (I+G)  <  II  G  |Lo_  (I  +  g"1) 

-  Ilg-ln  m  -  m  -2m 

4  * 

(b)  — i - o  (i  +  G)  <  a  (i  +  g'1)  <  || g_1  |L  a  (i  +  g) 

II G  ||2  ra  “  m  ~  2  m 


Proof :  Follows  immediately  from  Theorem  3  using 


I  +  G_1  E  G-1(I  +  G)  and  I  +  G  E  G(l  +  g"1) 
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For  the  rest  of  this  subsection  we  shall  let  H  denote  either 
I  +  G  or  I  +  G  ^  according  to  whether  additive  or  multiplicative 
perturbations  are  appropriate.  The  next  theorem  will  show  how  the 
singular  values  of  H  +  L  can  be  bounded  in  terms  of  |j  L  ||2  and  the 
singular  values  of  H. 


Theorem  5:  Suppose  0^(1!)  >_  >  0  for  some  k,  1  <_  k  m,  and 

||  L  8.  Suppose  further  that  B  <  a^.  Then: 

(a)  0.(1  +  H_1L)  >  1  -  f 

x  k 

(b)  (H  +  L)  _>  ak  -  8  . 

(Note:  If  k  ?  m,  H  +  L  is  not  necessarily  invertible  if  8  is  too  large.) 

Proof:  (a)  Use  1  =  1  +  H_1L  -  and  A  *  I  +  H*1!.,  B  =  -h'3'L,  i  =  k, 

j  =  m-k  +  1  in  Theorem  3(a)  to  get 

0  (I)  <  0  (I  +  H_1L)  +  0  ..  (H_:LL)  . 

m  —  k  m-k+1 

Thus  0.(1  +  H-1L)  >1-0  (H-1L) 

k  —  m-k+1 

>_  1  -  ||  L  jj^  *  Cm_k+1(H  1)  by  Theorem  3(b) 

-  1  -  II L  ll20k(H) 

>  1  - 

(b)  Use  H  =  H  +  L  -  L  and  A  =  H  +  L,  B  =  -L,  i  -  k,  j  *  1 
in  Theorem  3(a)  to  get 

Ok(H+L)  !<7k(H)  -  ||l||2  >  ak  -  8  . 
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The  case  k  3  m  is  of  special  interest  in  Theorem  5  as  it  bears 
directly  on  our  two  basic  inequalities  (18)  and  (21)  of  the  form 


which  are  sufficient  to  guarantee  stability  of  a  perturbed  closed-loop 
system.  Specifically,  if  ||  H  ||  2  <_  —  and  |( L  ||  <_  B  with  0  <_  &  <  a, 
then  H+L  is  invertible  and  ||(H+L)  ^||  <  ~—q  or  a  (H+L)  >  a  -  0. 

Note  that  Theorem  5  was  expressed  in  terms  of  isolating  |[  L  ||  -  By 
analogy  with  the  inequalities  (17)  and  (20)  we  can  also  have  the  fol¬ 
lowing  stronger,  but  perhaps  less  useful,  theorem. 

Theorem  6s  Suppose  a  ,  (H  ^L)  <1-6  where  0  <  6  <  1  and  1  <  k  <  m. 
Then: 

(a)  Ok(l  +  H_1L)  >  6 

(b)  a.  (H  +  L)  >  — ^r- 

k  ’  IIh-1|!2 

Proof :  (a)  From  the  proof  of  Theorem  5  we  have 

a.  (I  +  H-1L)  >  1  -  a  ,  (H-1L)  >  6 

(b)  From  I  +  H  =  H  1 (H  +  L)  and  Theorem  3 (b)  we  have 

ak(i  +  h-1l)  <  ak(H  +  l)  •  ||  H _1  !|  2 

whence  a,  (H+L)  >  - =■ -  . 

k  "  llH"1 1|2 
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d.  SPECIAL  RESULTS  WHEN  G(s)  =  C(sl  -  Aj^B 

In  this  subsection  we  shall  make  use  of  the  fact  that  the  frequency 
response  matrix  is  of  the  form 

G(j<o)  -  CtjMI  -  A)_1B 

Let  us  further  define 

F(ja))  =  C(j«l  -  A  +  BC)-1B  (30) 

Recall  the  Sherman-Morrison-Woodbury  formula: 

(W  +  XYZ)'1  =  W"1  -  W-1X(Y_1  +  ZW'1X)"1ZW~1 


assuming  the  indicated  inverses  exist.  Then  it  is  easy  to  verify  that 


(I  +  G(jai))'1  =  I  -  P ( j(i>) 


(31) 


and,  from  (27) , 


(I  +  G"1(jw))“1  =  F(jW) 


(32) 


Thus  our  results  in  the  last  section  (for  example,  Theorems  4,  5, 
and  6)  can  all  be  cast  in  terms  of  F  by  noting  that 


V1  +  G)  *  5 


m-k+1 


(I  -  F) 


(33) 


and 


V1*5'1’  -  -  .TTf) 


(34) 


m-k+1 


Moreover, 


||(I  +  «r1||  -  ||i-  F 


(35) 
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and 

)|{I+  g"1)"1!!  -  || F  ||  (36) 

for  any  of  the  norms  we  have  been  considering  (in  particular.  It  *  n  in 
(33)  and  (34)).  Use  of  (31)  and  (32)  results  in  an  apparent  savings 
in  the  number  of  linear  systems  to  be  solved  (i.e.,  number  of  inversions) 
and  we  shall  exploit  this  fact  in  the  next  section. 
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4.  COMPUTATIONAL  PROBLEMS 

a.  COMPUTATION  OF  FREQUENCY  RESPONSE  MATRICES 

As  we  have  seen  above,  an  object  of  considerable  interest  in  studying 

the  robustness  of  stability  of  linear  systems  is  a  graph  of  - — - 

||(I  4.  G"1(jco))‘1|| 

as  a  function  of  co.  When  G(jai)  =  C(j6Jl  -  A)  1B  we  saw  that  ||(I  +  G(j£u))  ^||  » 
||l  -  F  ( jto )  ||  and  ||(I  + G_1(ja)))”1|j  =  \\f  ( jco)  |j  where  F(jw)  =  C(jo)I  -  A  +  BC)-1B. 
Thus,  regardless  of  the  norm  used,  a  quantity  of  the  form 

C(jO)I  -  H)_1B  (37) 


must  first  be  computed.  We  shall  assume  throughout  this  and  the  next  sub¬ 
section  that:  (i)  B  e  3RnXm,  C  effi”11®,  H  e®™  are  given 
(ii)  n  ^  m 

(iii)  (37)  is  to  be  evaluated  for  a  large  number,  N,  of 
values  of  0);  typically  N  >>  n. 

Bather  than  concentrate  on  exact  operation  counts,  which  may  be  fairly 
meaningless  anyway,  we  shall  give  only  order-of-magnitude  estimates. 

It  will  be  seen  that  the  bulk  of  the  computational  load  rests  on  evalu¬ 
ating  matrices  of  the  form  (37)  and  so  we  shall  focus  initially  on 
that  problem. 

nxn  _i 

If  A  eiR  is  dense,  the  most  efficient  evaluation  of  C(jcoI  -  A)  B 

by  an  LU  factorization  of  A,  solution  of  m  triangular  systems  to  get 

(jcol  -  A)  ^B,  and  finally  a  matrix  multiplication,  requires  approximately 

13  12  2 

— n  +  —  ran  +  m  n  multiplications  (and  a  like  number  of  additions;  we 
shall  henceforth  count  only  multiplications).  This  figure,  when  multiplied 


by  N,  represents  a  rather  large  amount  of  computation. 


If  A  is  initially  transformed,  however,  the  computational  burden 
can  be  reduced  quite  considerably.  If  T  is  a  similarity  transformation 
on  A  we  have 

CtjCOl  -  A)_1B  =  CT  ( jcol  -  T  1AT)”1T_1B  . 

Let  us  define 

H  -  T_1AT 

and  agree,  for  convenience  to  still  label  CT,  T  1B  the  transformed  C  and 
B  matrices,  respectively,  as  C,  B  respectively.  We  now  have  the  problem 
of  evaluating 

C( jO>I  -  H)_1B 

where  H  may  now  be  in  such  a  form  that  (jcol-  H)-1  can  be  computed  in 

less  than  0(n3)  operations.  For  example,  A  can  always  be  reduced  to 

upper  Hessenberg  form  by  (stabilized)  elementary  transformations  (f-  n3 

6 

5  3 

multiplications)  or  by  orthogonal  transformations  (y  n  multiplications) . 

These  transformations  are  very  stable  numerically  and,  while  0(n3),  are 

performed  only  once  at  the  beginning  of  the  calculations.  The  resulting 

linear  system  to  be  solved  -  for  N  different  values  of  w  —  now  has  an 

upper  Hessenberg  coefficient  matrix  and  can  be  solved  in  approximately 
1  2 

—  mn  multiplications.  Moreover,  Hessenberg  systems  can  be  solved  very 

accurately  with  the  growth  factor  in  Gaussian  elimination  bounded  above 

by  n?  see  [8].  Computing  C(jo>I  -  H)  ^B  still  requires  an  additional 
2 

m  n  multiplications.  Neglecting  the  initial  transformation  and  deter¬ 
mination  of  CT  and  T  ^B,  the  Hessenberg  method  requires  approximately 
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12  2 

—  mn  +  m  n  multiplications  (for  each  value  of  to),  a  considerable  savings 
over  the  O(n^)  algorithm  if  n  »  m. 

Of  course,  other  transformations  T  are  possible.  One  possibility  is 
to  reduce  A  to  upper  triangular  (Schur)  form  by  means  of  orthogonal  simi¬ 
larities.  This  is  considerably  more  expensive  than  reduction  to  upper 
Hessenberg  but,  again,  need  only  be  done  once  at  the  beginning.  How¬ 
ever,  the  resulting  linear  system  to  be  solved  at  each  step  is  upper 

2 

triangular  and  so  still  requires  0(mn  )  multiplications.  Because  of 
potential  difficulties  with  multiple  eigenvalues  of  A  there  seems  to  be 
little  real  advantage  gained  by  this  procedure.  Substantial  savings 
could  be  gained  though  if  the  eigenstructure  of  A  were  such  that  it 
was  diagonalizable  by  a  reliably  computable  T.  Since  this  involves 

9 

consideration  of  the  essentially  open  numerical  problems  associated  with 

computing  invariant  subspaces  we  shall  not  pursue  the  details  here. 

But  assuming  such  a  transformation  were  possible,  C { jtol  -  D)  *3  with 

2 

D  diagonal,  could  be  computed  with  approximately  mn  +  m  n  multiplications 
for  each  value  of  to.  Attractive  as  this  appears,  the  potential  for  severe 
ill-conditioning  of  the  eigenproblem  associated  with  A  render  this  latter 
method  unreliable  as  a  general-purpose  approach.  We  shall  subsequently 
consider  only  the  Hessenberg  method. 

The  analysis  above  has  been  done  under  the  assumption  that  complex 
arithmetic  was  performed.  We  now  outline  how  G  ■  C(jtoi  -  H)  ^ B  might 
be  determined  using  only  real  arithmetic.  The  matrix  H  is  assumed  to 
be  in  upper  Hessenberg  form.  We  wish  to  solve  first 


(jtol  -  H)Z  -  B 


(38) 
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G  -  CZ 


Suppose  Z  ■  X  +  jY  where  X,  Y  eE^’lXin.  Upon  equating  real  and  imaginary 
parts  in  (38)  we  get  the  following  order  2n  real  system  to  determine 


X  and  Y: 


C :)  C)  ■  0 


1  22-1  22 
Thus  X  ■  —  HY  and  Y  *  — co (o>  I  +  H  )  B.  The  matrix  (to  I  +  H  )  will  be 
to 

2  2 

invertible  if  (jcol  -  H)  is  invertible.  Note  that  (to  I  +  H  )  is  no  longer 
upper  Hessenberg  but  is  almost  in  the  sense  of  having  two  rather  than  one 
nonzero  subdi agonal.  Its  shape  is  wholly  typified  for  n  =  5  by  the 
matrix 


X  X  X  X  X 


X  X  X  X  X 


X  X  X  X  X 


0  x  x  x  x 


ip  0  X  X  x 


Linear  systems  involving  matrices  of  this  type  can  be  solved  using 
2 

approximately  n  multiplications.  We  summarize  the  Hessenberg  method 
using  real  arithmetic: 

(i)  Reduce  A  to  upper  Hessenberg  form  H,  transform  B  and  C, 

2 

and  compute  H  ;  this  step  is  done  only  once. 


2 

(ii)  Solve  (to  I  +  v  {  =  -toB  for  Y. 

(iii)  Compute  X  *  ^  HY  . 

(iv)  Compute  G  ■  (CX)  +  j  (CY)  . 

2 

Step  (ii)  requires  approximately  mn  multiplications,  step  (iii)  requires 

12  2 
approximately  —  mn  ,  and  step  (iii)  approximately  m  n.  The  total  number 

3  2  2 

of  multiplications  is  approximately  —  mn  +  m  n. 

Storage  requirements  for  the  Hessenberg  method  with  real  arithmetic 
are  approximately  double  those  for  complex  arithmetic. 

b.  COMPUTATION  OF  ROBUSTNESS  MEASURES 

We  have  seen  above  that  quantities  of  the  form  (37)  can  be  reliably 
2 

evaluated  in  0(mn  )  operations.  There  then  remains  the  problem  of 
determining  (35)  or  (36). 

C««e  1;  ||  •  ||2 

For  (35),  the  singular  value  decomposition  (SVD)  of  I  +  G(jto)  can 
be  computed  for  each  value  of  to.  Each  SVD  typically  requires  approximately 
6m3  multiplications.  The  smallest  singular  value  is  then  the  quantity  of 
interest.  For  (36) ,  inversion  of  G  can  be  avoided  by  finding  the  SVD 
of  F ( jto) ,  again  in  approximately  6m3  multiplications .  The  inverse  of 
the  largest  singular  value  of  F  is  then  the  quantity  of  interest. 

Caae  2:  ||  •  1^  or  ||  •  ||  ^ 

Use  of  either  of  these  norms  in  (35)  or  (36)  involves  negligible 

2 

computation  as  compared  to  Case  1,  namely  about  m  additions  and  absolute 


values  and  m-1  arithmetic  comparisons . 


In  both  cases,  the  additional  work  required  is  usually  small  compared 
2 

with  0{mn  )  especially  if  n  >>  m.  However,  if  m  is  large  relative  to 
n,  significant  savings  cam  be  realized  in  using  [|*  ||^or  ||*  rather  than 
||  *  In  fact,  using  our  previous  approximate  operation  counts  for 

the  Hessenberg  method  and  setting  n  =  km,  we  have 


work  per  value  of  u  using  ||  *  |  \2  k2  +  2k  +  12 

p  «  "  ”■*  2 

work  per  value  of  w  using  ||  •  ||^  or  ||  *  1^  k  +  2k 

~  k^  +  2k  +  24 

Note  though  that  p  ~  - = -  if  singular  directions  are  also  com- 

x  +  2k 

puted. 

In  the  event  A  (or  A  -  BC)  can  be  successfully  diagonalized  as 
mentioned  in  Section  4. a.  the  potential  savings  in  avoiding  ||  •  ||2  are 
somewhat  greater.  In  fact,  we  then  have 


(or  p  *  — —  if  singular  directions  are  also  computed)  . 

The  above  comparisons  are  only  approximate  and  should  in  no  way 
be  construed  as  definitive  statements.  The  purpose  of  this  section  is 
to  merely  introduce  certain  aspects  of  the  numerical  computations  and 
suggest  further  avenues  of  exploration.  A  great  deal  of  numerical 
experimentation  remains  to  be  done.  Reliable  software  such  as  LINPACK 
[9]  for  linear  systems  will  be  of  great  benefit  in  this  research. 
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5.  CONCLUSIONS 

We  began  this  paper  with  an  attempt  at  a  "formal"  definition  of 

robustness.  We  then  applied  the  definition  to  the  problem  of  robustness 

of  stability  of  linear  systems  as  discussed  in  [1] .  The  cases  of  both 

additive  and  multiplicative  perturbations  were  discussed  and  a  number  of 

relationships  between  the  two  cases  were  given.  Finally,  a  number  of 

computational  aspects  of  the  theory  were  discussed  including  a  proposed 

new  method  for  evaluating  general  transfer  or  frequency  response  matrices. 

The  new  method  is  numerically  stable  and  efficient,  requiring  only 
2 

0(mn  )  operations  to  update  for  new  values  of  the  frequency  parameter 

3 

rather  them  0(n  ) 

A  number  of  interesting  research  areas  suggest  themselves  in  this 
work.  One  such  area  is  that  of  constrained  perturbations.  For  example, 
in  our  basic  problem  we  were  concerned  with  the  nearness  to  singularity 
of  a  nonsingular  matrix  A  €  E030"1 .  If  the  admissible  perturbations  E 
are  somehow  constrained  for  one  reason  or  another,  for  example  E  upper 
triangular,  the  usual  bound  on  ||  E  ||  for  which  A  +  E  is  singular  but  E 
is  "dense"  may  be  overly  pessimistic.  Related  to  this  is  the  fact  that 
our  bounds  were  derived  for  the  "worst  case".  The  size  of  perturbations 
allowed  in  a  linear  system  to  ensure  continued  closed- loop  stability  may 
very  well  be  larger  than  we  have  derived  if  inputs  to  the  system  are 
constrained  in  certain  directions. 

We  have  concentrated  in  this  paper  on  the  analysis  of  linear  control 

systems.  There  are  many  interesting  —  and  difficult  —  synthesis  problems, 

however.  For  example,  can  A,  B,  C  be  chosen  to  assign  certain  singular 
+1 

values  of  1  +  G'  ?  What  is  the  effect  of  changes  in  B  or  C  on  the 


njewBSa 
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±1  ±1 
behavior  of  .1  +  G_  ?  Can  a  matrix  K  be  determined}  so  that  I  +  (GK ) 

has  certain  singular  values? 

On  the  computational  side,  more  research  needs  to  be  done  on  updating 
parametric  problems.  That  is,  suppose  we  have  a  matrix  (say,  G ( jco) ) 
which  depends  "in  a  rank  m  way"  on  a  parameter  co.  When  oj  changes  how 
can  various  quantities  be  updated  efficiently? 

Finally,  as  mentioned  in  Section  4.b.,  a  great  deal  of  numerical 
experimentation  is  necessary  to  get  a  qualitative  feel  for  the  numbers 
in  determining  robustness  measures. 
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